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ABSTRACT 


The general design of a system of subroutines for solving the initial value 
problem in ordinary differential equations is given. An attempt has been 
made to design these subroutines in such a way that they will be easy to 
use on easy problems and stiH be flexible enough to treat any type of 
initial value problem with a high degree of efficiency. Emphasis is on the 
use of these subroutines, rather than on the mathematical algorithms which 
at this time are not completely specified. Implementation of our design 
in FORTRAN TV suffers from deficiencies in the design of the multiple 
entry feature provided in some of the current FORTRAN IV compilers. 
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I. HT'TRODUCTIOjW 

This paper gives the form in which we plan to provide subroutines 
for the numerical solution of ordinary differential equations. Our 
goal is to provide a reliable algorithm, which incorporates features 
of use in a wide variety of applications, which is highly efficient 
on a wide class of problems, and which is easy to use on easy problems. 
This document has two primary purposes, 

1, We are soliciting comments from users and numerical analysts on 
the following; features not included which someone may find 
useful, features provided in a way more awkward than necessary, 
and the general approach. 

2. We want to call attention to the sad state of the multiple entry 
feature provided in the FORTRAN compilers of third generation 
computers. None that we know of can be trusted to work as it 
should with the kind of usage outlined below. (The problem and 
some possible solutions will be sent on request.) Past experience 
with our UWrVAC llo8 indicates that with a little care (i.e, code 
which takes into account quirks in the compiler) we can get the 
desired results. On the IBM 3^0 things are a little worse, since 
in this case the end user will have to use tricks in his coding 
to get the desired results, and the higher the degree of optimi- 
zation, the trickier he may have to be. The CDC 6000 series 
provides a multiple entry in name only, and thus the usage given 
below is not even permitted. Clearly it is poor algorithm design 
to be so computer dependent, and yet the other choices are so 
unsatisfactory that we take the approach below in the hope that 
the future will find the most widely used scientific programming 
language is one which has a well defined multiple entry feature. 
(Anyone for PL/1?) 

In addition we hope that there may be some ideas here of use to 
others who are engaged in the design of mathematical software. 

The design given here extends and modifies the features available 
in the integrator described in [1], The main extensions are the 
special provisions for stiff equations and for variational equations. 
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The overall design has resulted from experience with the integrator 
in [1] and all of the features included here are of use in some 
application. The features are given in all of their gory detail, 
since we believe the detail necessary to the evaluation of how a 
feature will v;ork in practice. The reader should, of course, skip 
over the details of features in which he has no interest. Results 
obtained with the integrator in [1] are given in [2], 

Central to our approach is the use of reverse communication , which 
means that whenever additional information is to be communicated, a 
subroutine returns control to the program which called it. This is 
opposed to the common practice in integration subroutines of passing 
the name of the subprogram which computes the derivatives through the 
calling sequence of the integration subroutine, and then having the 
integration subroutine call the subprogram of this name whenever a 
derivative evaluation is required. Those who have experience only 
with the latter procedure may not be impressed with the facts that 
reverse communication may mean a little less fussing with C0MM0W, 
that it clarifies the structure of a large program, and that on 
simple problems only one program need be written. We adopted the 
reverse communication concept because Charles Lawson had made it 
the policy for subroutines in the JPL subroutine library. After 
trying both approaches we prefer to use reverse communication, and we 
believe that most people who give a fair trial to both approaches 
will feel the same way. 

The following section lists the subroutines, gives their 
function, their entry points, the arguments passed through each 
entry point, and the reason for each entry point. This skeleton 
outline is meant only to serve as a guide to Sections III and IV 
which describe how to use these subroutines. Also see Figure i 
which outlines how the various subroutines communicate. 

The remainder of this introduction gives the class of problems 
which the subroutines can be used to solve and gives a very brief 
sketch of the formulas which are used to effect the solution. The 
details here need not be completely understood in order to comprehend 
most of what follows. 
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The primary problem is the initial- value problem 
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The formulas which are used require backward differences v ' 
only for i=A> 0 £ X s d, where the choice of X determines the 
characteristics of the formula. The larger the value of X the more 
accurate the formulas, but because of stability considerations values 
of X < d, are sometimes best. If X < d, the corrector equations 
given below require (in general) the solution of a nonlinear system 
of equations for • • ♦ » and in the 

formulas are constants which can be determined by requiring that the 
formulas be exact for as high degree polynomial as possible. 

The predictors are: 
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If i=6. formulas (1.5), (1.7), and (1.8) are not used. If ^=?d-l, 
formula (1.7) is not used. If jJ < d, in formula (1,6) is 

replaced by Finally if z' 'in eq, (1.2) is replaced witb 

zero, then the right hand side of formula (1,8) is replaced with 
zero. Note that in this last case one must have & < d. 
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II. THE SUBROUTINES 


A. The Core Integrator — DV^A 

The core integrator has the following jobs. 

1. Keeping track of the various ways the user can request 
output points, and returning control to the user at these 
points with information of interest to the user set to the 
proper values. 

2. Keeping the past history of the solution, which is stored 
in difference tables, up to date, 

3. Predicting the new value of the dependent variables on each 
new step. (Formulas (1.4) and (1.5).) 

4. Correcting the values of those dependent variables which 
are not being treated as stiff equations. (Formula (1.6), 
if jj=d. If ^^d, all corrections are carried out in DSTF. ) 

5. Estimating the error with the current stepsize, 

6. Estimating what the error would be if the stepsize were 
doubled, 

7. Determining if the user has requested more accuracy than 
the precision of the computed derivatives warrants, 

8. Selecting the order of integration formula to be used on 
each equation in the system of differential equations, 

(The order is selected independently for each equation in 
the system, ) 

9. At the option of the user, to provide output which enables 
one to see why the integrator is selecting the stepsize and 
the integration orders the way it is. 

10, Selecting the stepsize. (That is, deciding when to halve 
and when to double the mesh, ) 

11, Changing the stepsize, and the difference tables to correspond 
to the new stepsize. 

12, On user option, to select the value of jf in formulas (1.4) 

to (1.8) when the current value of A=d. The user must provide 
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tile initial value of (This option V7ill only be made 
available if we can find a simple reliable way of providing 
it. ) 

The SUBR0UTINE statement has the form 

SUBR0UTIKE DV0A(HEQ,T,Y,F,KI),EP, IFLAG,H,HMIMA,DELT,TFIHAL, 
MXSTEP, KSTEP , KEMAX, ElxlAX, KQ, YN, DT ) 

One calls DV0A to set up initial communication and to initialize 
certain variables internal to DV0A. Other entries and their 
function are given below, 

ENTRY DV0AR 

DV0AR is called to restart an integration, while maintaining a 
distance of CELT between the output points that one gets with 
IFLAG=3. 

ENTRY DV0Ai6(LEX,LSS,L00,LH0) 

DV0A0 is called to set up certain special returns, and optional 
output . 

ENTRY DV0AS(KS,ASTIF,MCHANG,STIFAC ,SGAM) 

DV0AS is called by the subroutine DSTF, to set up communication 
between DV0A and DSTF. The user who doesn't want to set up 
code for stiff equations can call DV0AS in order to find out 
if DV0A finds some of his equations to be stiff, 

ENTRY DV0AV(NVE,EPSV) 

DV0AV is called to set up DV0A to handle variational equations. 
ENTRY DV0AH(HMAXA,LHC) 

DV0AH is called to specify a maximum steps ize, or to specify when 
the steps ize is to be halved or doubled. 

ENTRY DV0AG(IG) 

DV0AG is called by DAGS to take care of the various disruptions 
that occur when locating GSTOPs, 

ENTRY DV0A1 

DV0A1 is called by the user whenever he wants to get back to the 
integrator and simply continue with the integration. 
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B. The Interpolator — DAINT 

The interpolator computes the values of the dependent variables 
and their derivatives (up to the order of the differential 
equation) at arbitrary values of the independent variable. 

The SUBR0UTIKE statement has the form 

SUBR^SUTIWE DAIHT(T,Y,F,KD,KQ,YNjIlT,WEV,IIY,TL,H) 

One calls DAINT to set up initial communication. This set up 
call is made by DV0A if one is doing an integration, 

ENTRY DAINTl 

DAINTl is called whenever one wants interpolated values to be 
computed. 

C. The GSTOP Locator — DAGS 

The GSTOP locator finds the zeros of arbitrary functions 

g-(t,y), where t and y are as defined in eq. (1.1). 

0 

The SUBR0UT1RE statement has the form 

SUBROUTINE DAGS ( T , Y , IFLAG, NG, NST^P, LTGS , G , GT ) 

One calls DAGS to set up initial communication, 

ENTRY DAGSl 

DAGSl is called to check for zero crossings, and also in the 
process of locating the zeros. 

D. The CHECK Subroutine — DCHK 

The CHECK subroutine gives a simple way for the user to get 
labeled output of everything in the calling sequence of DV0A, 
'Such output is frequently useful in debugging a program. 

The SUBROUTINE statement has the form 

SUBROUTINE DCHK(NEQ, T, Y , F , KD , EP, IFLAG, H, HMINA, DELT , TFIHAL , 
MXSTEP, KSTEP, KEMAX, EMAXjKQ, YN, DT, lO) 
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One calls DCHK to set up initial communication. 

ENTRY DCHKl 

DCHKl is called whenever one wants to print out the values of 
the variables in the calling sequence of DV0A, 

E, The Stiff Integrator — DSTF 

The stiff Integrator has the following jobs, 

1, Form the matrix used in the iteration for solving eqs, (1.6)- 
(1.8), This matrix is described in subsection. 17 E below, 

2, Deciding when to compute the Jacobian matrix. 

3, Carrying out the iterations used to solve eqs, (l,6)-(l.8). 

U, At the option of the user, to provide output of internal 
variables associated with the solution of eqs, (l.6)-(l,8). 

5. On user option, to select the value of A used in eqs. (1.4)- 
(1.8), when f < d. 

6. Changing the difference table from one value of A to the 
corresponding table for a different value of 1 . 

7. Obtaining an estimate of what accuracy can be achieved in 
solving eqs. (l.6)-(l,8). (This will only be done if it is 
-necessary for job 7 in the core integrator.) 

8. Determining if the user has computed the Jacobian incorrectly. 

The SUBROUTINE statement has the form 

SUBROUTINE DSTF(NEQ,T,Y,F,KD,EP,JFLAG,H,HMINA,DELT,TFINAL, 
MXSTEP , KSTEP , KEMAX, EMAX , KQ,, YH, DT , KB , AST IF , JSTOR , 
FPS,IFPS,WS,OTS,IW) 

One calls DSTF to set up initial communication with DSTF and DV^A, 
and to initialize certain internal variables in these subroutines. 
Other entries and their function are given below. 

ENTRY DSTFR 

DSTFR is used to restart an integration using DSTF in exactly the 
same way as DV^fAR is used with DVOA. 
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ENTRY DSTFj5(n®’,L0l) 

DSTE0 is called to set up a special return just before the matrix 
used in the iteration is factored, and to set up optional output. 

ENTRY DSTFC 

DSTFC is called when the user wants to change the values of jj 
use'd in eqs. (l.4)-(l,8), 

ENTRY DSTFJ 

DSTFJ is called if one wants to use a new Jacobian matrix in the 
iteration to solve eqs. (l.6)-(l.8) without the Jacobian calculation 
being requested by DSTF. 

ENTRY DSTFV(NVE,EPSV,IREF) 

DSTFV is called to set up DSTF to handle variational equations. 

The Partial Derivative Generator — DPART 

The partial derivative generator is used to generate the Jacobian 
matrix required by the stiff integrator. It is used when the 
user doesn't want to (or can't) write the code for the analytic 
partial derivatives. More work is required before specifications 
for this subroutine can be given. 
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III. USAGE 


A, The Basic Framework 

In this subsection we describe how to use the subroutines to 
integrate the basic initial- value problem (l.l). We begin by- 
listing all of the parameters in the call to the integrator 
together with how they are used. This is followed by a sample 
FORTRAN-like program which can be used as a guide when writing 
a program which calls the integrator. In the sample program we 
use symbols for statement numbers. Some of these symbols are 
referred to later in connection with other options, 

NEQ number of differential equations in the system, (Same 

as m in eq, (l.l). ) 

T independent variable, (Same as t in eq. (l.l),) 

Y vector of dependent variables. (Same as y in eq. (l.l).) 

For a system of first order equations Y(I) is the I-th 
dependent variable. In the case of higher order 
equations, Y(l) is the same as the I-th component of y 
as indicated just below eq. (l.l). Thus for the system 
F(l) = U", F(2) = V"; Y(l) = U, Y(2) = U% Y(3) = V, 

Y(4) = V. 

F vector of derivative values, (F(I) is the same as f^ 

in eq. (l.l),) The user must provide code to compute 
F given T and Y. 

KD order of the differential equations in the system. 

Thus for a system of first order equations, KD=1; for 
a system of second order equations, KD=2; etc. If 
differential equations of mixed order are to be inte- 
grated, KD must be a vector and KD(l) < 0 to inform 
the integrator that this is the case. The order of 
the I-th equation is then given by jKD(l)]. (Thus 
|KD(l)l=dj in eq. (l.l); for I > 1, KD(l) may be 
either positive or negative. ) Differential equations 
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with order greater than 4 can only be integrated 
by breakiag them into lower order equations, or by 
changing certain DATA statements in the subroutines, 

EP bound on the estimated local error. The local error 

in integrating each differential equation is kept less 
than EP/lO unless noise appears to be limiting the 
precision. If noise appears to be limiting the 
precision then the local error estimate is permitted to 
exceed EP/lO; for further details on this, see the usage 
of EMAX and IFLAG=6. For different error bounds on 
different equations, let EP(l) < 0 for I < K, and let 
EP(K) s 0. The local error control for the I-th 
equation is then based on |EP(l)j for I < K and on 
EP(K) for I s K. For differential equations of order 
d > 1, the error estimate is for the error in the (d-l)-st 
derivative. In this case the value of EP necessary to 
maintain a given accuracy in all lower order derivatives 
depends on the scaling. For a relative error bound, see 
subsection III B below. EP = 0 is not permitted except 
when HMAXA = 0. For further details on this see 
subsection III G below. 


IFLAG parameter used for communication between the integrator 
and the user. The integrator sets IFLAG as follows: 

=1 The value of Y for the current step has been predicted. 
The user should compute F and call DV0A1. 

=2 The value of Y for the current step has been corrected. 
The user should compute F and call DV0A1. 

=3 An output point has been reached (see the usage of DELT), 
to continue the integration call DV0A1. 

=4 T=TFIRAL, To continue the integration with a different 
TFUIAL change TFUSAL and call DV0A1, If DV0A1 is called 
without changing TFINAL, IFLAG is set equal to 8. 

=5 KSTEP S KS0UT. (See the description of MXSTEP. ) 
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=6 EMAX > ,1 and it appears that reducing H will not 

help reduce the global error because of round-off error 
or noise in computing F, A larger value of EP (or of 
|ep(KEMAX) j if EP is a vector) should probably be used. 

If EP is not increased, too small a steps ize is 
liable to be used. We have found that replacing EP with 
32*EMAX*EP works reasonably well. (Note that in most 
cases, EMAX will be only slightly larger than .1.) 
Increasing EP in this way should not degrade the accuracy; 
however, if the nature of the problem changes it may 
pay to use a smaller value of EP later in the integration. 

=7 |h| < HMINA. This may be the result of halving H, of 

coming to the end of the starting phase with jHj < HMINA, 
or of the user increasing the value of HMINA, If 
one wishes to continue with the current value of H, 
set HMINA ^ H and call DV0A1. If H has Just been halved 
(in which case EMAX > .1), one may continue with the 
old stepsize by simply calling DV0A1, (Such action is 
risky without a careful analysis of the situation. ) If 
the stepsize has not Just been halved, calling DV0A1 
will result in the integration continuing with the 
current value of H, and a return to the user with IFLAG=7 
will occur at the end of every step until )h| s HMINA, 

=8 Fatal error. The value of KEMAX gives the source of 
the error as indicated below. If DV0A1 is called an 
error message which includes the value of KEMAX will be 
printed, and the program stopped. 


KEMAX =1 KEQ ^ 0, 

=2 H=0 

=3 KD (or some |kd(I)|) is equal to 0 or greater than k. 
=4 BELT 0 (After IFLAG=3, Old T0UT=(Old T0UT)+DELT.) 
=5 TFINAL was not changed after IFLAG=4. 

=6 EP (or some EP(I))=0. 
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~7 At the initial point the direction of integration is in 
the direction opposite from, TFINAL, and the estimated 
error in extrapolating hack to TFINAL from the initial 
point is too large, 

=8 During the integration, either TFINAL or DELT has been 
specified in such a way that the next output point, 

TMARK, satisfies (TL-TMAEK)/H > 8 where TL is the last' 
value of T that the integration has reached. (Thus 
the output point specified is one that the integration 
process passed over 8 steps ago, assuming that H hasn't 
changed.) If one wants to continue under such conditions, 
set IFLAG=0, call DV0A1, and Y and F will be computed 
at this output point and a return will be made with 
IFLAG=3 or 4 as usual. Otherwise a call to DV0A1 results 
in the program being stopped as indicated above. 

=9 DV0A1 was called immediately after a call to DV0AV or 
DV0AS. After a call to DV0AV or DV0AS, either DV0A 
or DV0AR must be called. 

=10 NVE is specified improperly (possible only when the 
entry for variational equations is used), 

— =11 DV0AH was called with LHC0O and IFLAG09. 

=12 After a return from DV0A, DAGSl was called with IFLAG 
equal to something other than 1, 2, 4, or 9. 

=13 DAGSl was called after DAGS was called with NG ^ 0. 


The following errors are possible only with stiff equations, 

KEMAX=14 M0D( |ks(i) ] ,10). > 4 or KS(I) = -k*lO where k is a 
positive integer. 

=15 JST0R specified improperly. 

=l6 IFPS is too small. 

=17 INS is too small, 

=l8 DSTFl was called immediately after DV0A set IFLAG equal 
to something other than 1, 2, or 10. 
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HMINA 


DELT 


■19 DSTFl was called immediately after a call to BSTFV, 

After a call to DSTFV either DSTF or DSTFR must be 
called. 

■20 IREF (in DSTFV entry) is specified improperly. 

:21 The value of KS for a variational equation does not 

agree with the value of Iffl for the equation with which 
the variational equation is associated 

the stepsize. The initial value of H determines the 
direction of integration. The following points should 
be considered in selecting the initial value of H. 

1. Later values of H=2 •(initial H) where k is an 
integer. 

2. Efficiency does not depend critically on the initial 
value of H, see for example Table 12 in [2]. 

3. If it does not lead to problems in computing the 
derivatives (e.g, because of overflow or because of 
trying to compute the square root of a negative 
number), it is better to choose H much too large 
than too small. 

4. If one wants to obtain output values without the 

integrator doing an interpolation, then one must choose 

Ic 

H and DELT such that DELT=H*2 , k a non-negative integer. 

(In addition DELT must be an integer times a power of 

2 since otherwise there is liable to be a problem 

-5 -7 

with round-off error. Thus DELT=2 3x2 13x2 , 

-Q 

51x10 are all satisfactory, but DELT=.l will 
always lead to a need to perform interpolations.) 

minimum stepsize permitted after the integration is 
started. The user cannot restrict the minimum value of 
]Hj while the integration is getting started. After 
the integration is started, and whenever H is halved, 

|h| is compared with HMIM. li* Ih| < HMINA control is 
returned to the user with IFLAG=7. 

output increment. Initially the integrator sets the 
internal variable T0UT equal to the initial value of T. 
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TFIWAL 


MXSTEP 


KSTEP 

KEMAX 


EMAX 


Whenever T=T0UT a return is made to the user with IFLAG=3. 
If DV0A1 is called after IFLAG=3, T0UT is replaced with 
T+DELT, If T0UT does not fall on an integration step, 
interpolated values for Y and F are computed on the 
first step that (T-T0UT)*H > 0 and T is set equal to T0UT. 
If one selects an initial value of PELT such that 
H*DELT <0, the sign of PELT is changed hy the inte- 
grator and the initial value of T0UT is set equal to 
T+PELT. (Thus if H*PELT < 0 initially, one does not 
get output at the initial point.) 

final value of T. Vhien T reaches TFINAL, control is 
returned to the user with IFLAG=i|. Output values at 
TFINAL are obtained in the same way as described for 
T0UT in the description of PELT above. 

maximum 'number of steps. Initially the integrator sets 
KBTEP=0 and the Internal variable KS0UT = |MXSTEPj. 

KSTEP is increased by 1 at the end of each step. 

Whenever KSTEP S; KS0UT, a return to the user is made 
with IFLAG==5. If PV0A1 is called after IFLAG is set 
equal to 5, KS0UT is replaced by KSTEP+ | MXSTEP ] . If 
MXSTEP < 0, KS0UT is also replaced by KSTEP+ [MXSTEP } 
after IFLAG is set equal to 3 or 4, Thus if MXSTEP < 0, 
jMXSTEPj is the maximum number of steps that can occur 
before IFLAG = 3> 4 or 5. 

number of integration steps taken. KSTEP is initialized 
and incremented by the integrator, 

index of the equation responsible for EMAX (see below), 
KEMAX also indicates the type of error when IFLAG=8 
(see above). 

largest value in any equation of (estimated error )/e 
where g = jEP| for the equation under consideration. 
Ordinarily the stepsize is halved if EMAX > .1. However, 
with a recent history of round-off error or noise in 
the derivative evaluations limiting the precision, values 
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of EMAX as large as 1 are permitted. If round-off 
or noise appears to be limiting the precision on the 
current step, the stepsize is not halved under any 
c ir cums t anc e s . 

KQ, vector used to store integration orders. KQ(l) gives 

the value of q for the I-th equation. See equations 
(1.4) and (l»6) for the definition of q(i?=d). 

YN vector used to store Y at the end of each integration • 

step, 

DT a two dimensional array used to store the difference 

tables. DT(I,J)=V^“^^(J) where V^F^( j)=F^( j)= value 
of the J-th component of F on the n-th step, and 

V^'^\(J)=V^n^'^^‘’^n-l^'^^* k=l,2,...,KQ(j). 

The sample program is given below. 

In the type and dimension information given below, (?) indicates 
that the variable may be either a scalar or vector (see description 
above), m s NEQ, and n s total order of the system, 

.INTEGER NEQ,KD(?),IFLAG,MXSTEP,KSTEP,KEMAX,KQ(m) 

REAL EP(?),HMIWA,EMAX 

DOUBLE PRECISION F(m),DT(20,m) 

DOUBLE PRECISION T,Y(n) ,H,DELT,TFINAL,YN(n) 

Aq Assign the values for NEQ,KD,MXSTEP,EP,HMIM,DELT, and 

TPIKAL. Assign initial values for T, Y, and H. (These 
will be updated by the integrator. ) Values need not be 
assigned for IFLAG,KSTEP,KEMAX,KQ,EMAX,F,DT, and YN. 

A^ If certain special options are used they should be set 

up at this point, 

CALL DV0A ( NEQ, T , Y , F , KD , EP , IFLAG, H, HMINA , DELT, TFINAL , 
MXSTEP, KSTEP,KEMAX, EMAX,KQ, YN,DT) 

G0 T0 A^ 
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CALL DV0A1 

A^ G0 T0 IFLAG 

C0NTINUE (Special code is required here if certain 
options are used. ) 

Compute F(I), 1=1,2, ... ,NEQ (F=F(T,Y)) 

G0 T0 A^ 

Print results (T=T0UT=Last T0UT+DELT. See description 
of BELT.) Change CELT if desired (or ST0P if jDELTj 
was set large) 

G0 T0 A^ 

Print results (T=TFINAL) 

ST0P or Change TFINAL and 
G0 T0 A^ 

1^ Print results (KSTEP s KS0UT See the description of 

MXSTEP. ) 

G0 T0 Ag or ST0P (Depending on how you wish to use 
MXSTEP. ) 

Increase EP as s\:iggested under IFLAG=6 above. 

Print new EP, T, etc, 

G0 T0 A^ 

Print "H IS T00 SMALL" 

ST0P 

Ig G0 T0 A^ (There IFLAG and KEMAX will be printed and the 

program stopped. ) 

Note that the actions taken in the above program from 
statement to the end are but suggestions of what one might 
do for the specified values of IFLAG. Many other things can 
be done. 

The above basic framework may be added to as indicated in 
the options which follow. 
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B. Use of a Relative Error Test 

Users frequently desire to bound the local error by le*E(T,Y)| 
where e is a constant and E is some specified function of T and Y, 
This option is obtained in the program given in subsection A above 
with the following changes. At A^ assign a value to e instead of 
to EP; at I, set EP=e»E(T,Y) before continuing to and at 
increase the value of e instead of the value of EP, Either e or E 
can be a vector in which case EP is a vector, and the values stored 
in EP must be negative as mentioned in the description of EP 
previously. Users frequently set E equal to the Vector (|y(1)|, 
1y(2) 1 , . . . , Iy(NEQ) j ) (when KD=l) in order that the local error 
requested be relative to the solution. Unfortunately, this is 
frequently a poor choice. 

C. Restarting an Integration 

Because of discontinuities, or a change in the differential 
equations being integrated, it is sometimes best to restart the 
integration. This of course can be done by making the appropriate 
adjustments, calling DV0A and then going to A^. This procedure, 
however, is liable to change the spacing one gets on the output 
when IFLAG=3 because of the way T0UT is initialized, as explained 
in the description of DEBT, To preserve the spacing simply 
substitute 

CALL DV0AR 

for the CALL DV0A(IHEQ, . . , and then go to A^. 

D. Special Returns and Optional Output 

At any place in the program one may insert 
CALL DV0A0(LEX,LSS,L00,LH0) 

The variables LEX, LSS, L00, LH0 are all type INTEGER; and each 
refers to a special option as indicated below. Any option one 
is not interested in will be left unchanged if the corresponding 
variable is equal to 0. The integrator does not change the value 
of any of these variables, and any time the user changes one of 
them DV0A0 (or DV0A or DV0AR) must be called for the change to 
take effect. 
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LEX determines whether output at TFUJAL is to be obtained 

using extrapolation or interpolation. Extrapolation is 
useful if F has a singularity at TFIML, 

>0 Obtain values with extrapolation, (if TFINAl falls on 
an integration step, the extrapolation occurs before the 
step is taken. ) 

=0 Wo change. (The nominal internal flag is set for 
interpolation. ) 

<0 Obtain values with interpolation. 

LSS is used to specify whether one wants the integrator to 

return control at the end of each step. Control is 
returned with IFLAG=9. Thus if this option is used one 
must add an in statement and then add a statement 
mmbered where one does whatever he wants to do at the 

end of each step, followed by a 00 T0 A^. 

>0 Return control to the user with IFLAG=9 at the end of 
every step. 

=0 Wo change. (The nominal internal flag is set for no 
return. ) 

<0 Wo output. 

L00 is used to specify that the integrator is to print certain 

intermediate re stilts on every step. This information is 
useful for debugging and in analyzing the error control 
that is used. The output to be provided is not yet 
completely specified. 

>0 Give the output. 

=0 Wo change. (The nominal internal flag is set for no output.) 

<0 Wo output , 

LH0 is used to specify whether one wants control returned 

when the stepsize is changed. Control is returned with 
IFLAG=11 just before the stepsize is halved and with 
IFLAG=:12 just before the stepsize is doubled. Thus if 
this option is used one must provide for larger values of 
IFLAG in statement A^, Wo special returns are 
made if the stepsize is reduced in the process of 
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taking the first step. Control is returned with the 
current value of Y contained in YN and the current 
value of F(I) in DT(1,I). When the stepsize is halved 
Y and F contain the values that resulted in the decision 
to halve the step. To continue, call DV0A1, 

>2 Control is returned just before the stepsize is changed. 

=2 Control is returned just before the stepsize is doubled. 

=1 Control is returned just before the stepsize is halved. 

=0 No change. (The nominal internal flag is set for no 

return when H is changed.) 

<0 Control is not returned when the stepsize is changed. 


E. Variational Equations (and Others With Similar Characteristics) 

Consider the system of differential equations (compare with 
eq. (1.1)) 

Zi = i=l,2,...,m 

(■JJ 

Zi = i=m+l, . . . ,m+y, 


where all variables common to eq. (1.1) are defined the same as 


m+1 


m+2 


3 » • • > 


there, =. (?,^^) where •’ 

z ^ ) . Thus y. is formed by augmenting the vector y with y , 

A s 

the variables introduced in the bottom equation of (3E.1). If 

there are subexpressions in f^, i > m, which depend only on t 
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and y then the option described here permits one to compute 
these subexpressions only once per step instead of the twice 
per step that would ordinarily be required. In addition accuracy 
is slightly better since only corrected values of y are used in 
computing the subexpressions. A common situation where this 
option can save significant time is in the integration of 
variational equations, in which case the bottom equation in (3.1) 
has the form 


= A(t,y)y + b(t,y) 


(3E.2) 


where A(t,y) is a block diagonal 

matrix of the form 


'o'(t,y) 


0 


A(t,y) = 


a(ljy) 




Q-Ctjy) I 


> a(t,y) = 


m 


I 

« 

. , 'h 

« 

• 

dyn 

• 

• 

# 

df 

m . 

» 

df 

. . m 

5^1 

dyn 

1 

: d.. 

a ^ 

Also d . 
d 


, (3E.3) 


where y. is the i-th component of y and n = Z d.. Also d. = d. 

^ ^ t] ^ 


if jsi mod' m. One might have for example “ 

l,2,,..,n, where T] is the initial condition for y, , In this 
k -+ k 

case the first mn components of b are 0. If f.(t,y) can be 
written as f^(t,y,g^, . , . ) where the p's are parameters which 
are fixed on a given integration and one wants to compute 
5yj^/9Pj^ in addition to the Sy./r'Ilj^ computed as indicated above, 
simply sat = Syi/SS^ and - af 
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where b is the v~th component ov b. The matrix A then has 

V 

m(n+j) rows and n(n+j) columns, with the block diagonal structure 
as indicated in ( 3 E, 3 ). lu the case of variational equations 
this option cuts in half the number of evaluations of ^(tjy) and 
a considerable savings in some cases. 

This option requires the definition of two new parameters. 

KVE is the number of additional equations. (WVE = (j, and 

NEQ = m in (3.1) above. The total number of equations 
is given by lUEQ+NVE. ) 

EPSV is an error tolerance for the variational equations. 

If EPSV > 0, EPSV is used for EP in all equations with 
index > KEQ. If EPSV = 0 no error estimates are computed 
for those equations with index > NEQ, and the step- 
size will be selected on the basis of the estimated 
errors in the first KEQ, equations. If EPSV < 0, then 
|ep(I)|, (I > IffiQ) will be used in the error control 
on the I-th equation where, as before, if EP(l) s 0 for 
some I, this value will be used for all larger I. (Note 
that if EPSV > 0 and one is using a relative error 
test as described in subsection B, then EPSV must be 
changed in addition to EP when IFLAG=1. ) 

As implied in the above, Y is redefined now as the augmented 
vector y. in eq. ( 3 E.I) and F is extended to contain f., 

■A. IL 

i=l,2,...,NEQ+WVE. 

To use this option the following changes should be made in 
the sample program given in subsection A above. 

In the dimension information one must have m s NEQ+NVE, and 
n s total order of the system including the equations which have 
been added. If one has mixed orders, the dimension of KD must 
be s m, otherwise KD can be a scalar. One should add NVE to the 
variables specified to be of type INTEGER, and add EPSV to those 
specified to be of type REAL. 
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At Aq initial values for WVE and EPSV must be assigned. 

At A^ (before the call to DV0A) 

CALL DV0AV(UVE,EPSV) 

Code from A^ to the G0 T0 A^ just below should be replaced 
with; 

A^ G0 T0 IFLAG 

C0NTINUE (Special code is required here if certain 
options are used.) 

Compute F( I) , 1=1,2, , . . ,WEQ 

G0 T0 Ag 

I,* Compute subexpressions for extra equations. (o'(t,y) 

and b(t,y) in. the case of variational equations.) 

Compute F(I), I=NEQ+1, KEQ+2, ..., EEQ+KVE 
(Note that if the F(l) computed at depend only on 
y (thus not on y ) in eq. 3E.1, then these F(l) need 
not.be computed when IFLAG=2, since the values will be 
the same as those obtained when IFLAG=10. ) 

IF(IFLAG-2) Einy, V^, A^ ( IFLAG < 2 should be impossible.) 

If one is not using the end of step return (lSS) in subsection D, 
simply set Iq=A . 

y 2 

We have Introduced the feature described in this section in 
its simplest form for convenience in its description. It is 
actually available for the more general case 

(d.) 

^i ^ = fj^(t,y) i=l,2,...,m, 

i=m+l,...,m+]o,^, 

(d.) ^ (3E.4) 

= f^(t,y^) i=m+p,j^+l,. . ,,m+|j,j^+p,2, 

^ i " ^ ^ i=m+ix^+p,2+l , . . . , , 
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where = (y^y^) as before, = (y^jyi^) where y^ gives the 

the variables introduced in the third set of equations, y^ = 
(y-ciVn)) a'tc. To get this extended capability, KVE and EPSV 
are used as follows. 


NVE is a vector, with 


NVE(1) = -V where v is the number of additional sets of equations, 
and KVE(l) < 0 is necessary to indicate that NVE is a 
vector. Thus if one has the system given in eq. (3E.4) 
down to and including z, = f.(t,y ), then v=3 and 
NVE(1) = -3, and the total number of equations is 


m+(j,^+U2+U2 • 

NVE(2) is used to tell the user what set of equations requires 
a derivative evaluation when IFLAG=10. If we let = 
m=NEQ, then when IFLAG=10 and NVE(2)=j, F(l) is to be 
computed for I = |iQ+d'^+* * * *’*’l^j* 

NVE(2+k) is the number of equations in the k-th additional set 
(=^ji ), k=l,2,.,.,v. The total number of equations is 
thus given by NEQ+NVE(3)+" *+NVE(2+v) where v = -NVECI). 


EPSV is a vector, with 


EPSV(l) used to flag the usage of EPSV as follows 

>0 e . = EPSV{2) 

0 

=0 e. = 0 (In this case EPSV can be a scalar.) 
0 

<r 0 e. = EPSV(j+l) 
d 

Where 


> 0 



means use p. . for the error tolerance in every equation 
of the j-th additional set. (1=^,^+* • *+p,^_^+l, . .. .+ 

itj). 

means do not check the err.or in the j-th additional set. 
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s . < 0 
a 


means refer to EP(l) (passed through, the call to DV0A) 
for the error t oleranc e , j -1^^ j • • • ^ IJ-q"*" • * • ’’'M' j ♦ 

(of course, if EP(l) a 0 for some I in this range, that 
value will be used for all larger values of I in the 
j-th set of additional etiuations. ) 


Except for the obvious changes, usage is identical with that 
given earlier for the case when KVE is a scalar down to the 
statement labeled (see the middle of page 2h), From there on 

one should use the following. 


Compute subexpressions for the j-th additional set, 
where j=NVE(2). 

Compute derivatives for the j-th additional set. 

G0 T0 Ag. 

Compute F(I), for all I > HEQ. (As before, one should 
not bother to compute those F(l) whose values have not 
changed since they were computed when IFIAG=10.) 

G0 T0 V^. 


F, Direct Interpolation From Results, Present and Past 

During an integration one can obtain values' of Y and F at 
any time for any value of 'f by simply setting T to the desired 
value and then executing 

CALL DAINTl 

However, if one uses this option when IFLAG is not equal to 3? 
4, 5j 6, or 9 j then the values of T, Y and F must be saved and 
restored if one intends to continue the integration. Note also 
that one can expect reasonable accuracy in the interpolation 
only if 


-%.in ^ ^ ^ 


(3F.1) 
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where TL is the value of T at which the difference tables DT were 
computed, H is the stepsize, and = niih{KQ(l)}. 

To save the solution for retrieval at a later time the 
following changes should be made in the sample program of sub- 
section A, 

Just below Aq set an integer variable KST^RE = -1, at A^ insert 
CALL DV0A0(LEX,LSS,L00,LH0) 

where LSS > 0, LH0=1 (or is s 3 ) and subject to these constraints 
one has set the parameters passed in this call to appropriate 
values depending on what one wants. (See subsection III D. ) 

At A_ add (at least) I , L , and to the computed G0 T0. 

(If there are no variational equations, set I^q=A 2 .) At I^ the 
following code should be inserted. 

I^ Do any other thing you may want to do when IFLAG=9. 

KQMIN=100 

D0 1=1, NEQ, (Or HEQ+KVE if variational equations are 

present, ) 

KQMIN=MIKO( KQJ4IN, KQ( I) ) 

IF (KSTEP.LT.(kST 0RE+KQMHT)) G0 T0 Ag 

Ij^j^ Store T,H,YTT,KQ,DT (on tape, disc, or drum) 

KST0RE=KSTEP 

Make special test for IPLAG=11 and transfer if you want 
to distinguish between IFLAG=11 and IFLAG=9- 
G0 T0 A^ 

To retrieve the solution at a later time, one first assigns 
values to TL,H,YW,KQ,, and DT, from the values stored earlier, 
where TL is one of the stored values of T, and the remaining 
variables were stored at the same time as TL, The value of TL 
retrieved should satisfy (3F.1), where T is the value of the 
independent variable where values of Y arid F are desired. In 
addition to variables already defined, one must assign values to: 
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NEV v;hich is the total number of equations (=KEQ if there are 

no variational equations), and 

NY which is the total number of components in Y (=the total 

order of the system) 

Before the first interpolation, one makes the set-up call 

CALL daint(t,y,f,kd,kq,yu,dt,m;v,ny,tl,.h) 

(All of these variables have been defined above.) Then to obtain 
Y and F given T 

CALL DAINTl 

remembering that TL,H,YN,KQ, and DT are equal to previously stored 
values, and are selected as a function of T. 

One may not use DAIKT in this way at the same time as an 
integration using DV0A is being performed. One may of course, 
get such a capability by using another interpolation program 
identical to DAINT, except with the names DAINT and DAINTl changed 
to something else. (One may also have to change the names of 
other variables to avoid conflict with those used in the integration.) 

Of course, if one is only interested in interpolating for 
the variables associated with some of the equations> one need 
store only the corresponding values of YN,KQ, and DT. The 
obvious changes in KD,HEV, and NY shoiild then be made when 
retrieving the information. 


G. Getting More Direct Control Over the Integration 

The options described here give the user several levels of 
more direct control. One can merely specify that the steps ize is 
not to exceed a certain maximum value, or one can force the 
integrator to halve or double the steps ize at specified times, 
or one can specify the integration orders to be 
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used in addition to specifying the stepsize. These options are 
provided for those masochistic users who like to have enough 
rope to hang themselves. 

The following two parameters are used. 

HMAXA determines an internal variable, HMAX, so that doubling 
H is not permitted if doing so would result in an H 
which satisfies |h| > HMAXA. (Normally HMAX is so large 
that it never limits the stepsize.) HMA^ is type HEAL. 

LHC determines the action to be taken as follows. LHC is 

type INTEGER. 

=0 Assigns the internal variable HMAX as specified by 
HMAXA. 

=-l Same as LHC=0, except in addition the stepsize is halved. 

=-2 Same as LHC=0, except in addition the stepsize is doubled. 

=1 The stepsize is halved, and HMAX is set so that H will 
not be doubled unless the user so requests, (HMAXA is 
ignored. ) 

S2 The stepsize is doubled and HMAX is set so that H will 
not be doubled again unless the user so requests. 

(HMAXA is ignored. ) 

One gets the desired’^action -by means of the statement 
CALL DY0AH(HMAXA,LHC) 

This call can be made at any time if LHC=0. Otherwise it can 
only be made after IFLAG has been set equal to 9» Thus one must 
make use of the LSS option of subsection III D. If LHCj^O, one 
should call DV0AH only after everything else that one might want 
to do when IFLAG=9 has been done, and the call to DV0AH should 
be followed by a G0 T0 A^ in the sample program of subsection A. 

If LHC=0, one should simply continue as usual after the call to 
DV0AH, If DV0AH is called with LHC^O, and IFLAG was not previously 
set to 9} then after the return from the call to DV0AH, IFLAG will 
equal 8, a fatal error condition. In using the above options, one 
can insure that the stepsize will not be halved by setting EP 
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(or EP(1)) to a ridiculously large positive value, or by setting 
HMIWA to an appropriate value whenever the stepsize is changed. 

One can eliminate any error estimation, stepsize control, or 
integration order selection on the part of the integrator by 
calling DV0AH with HIyIAXA=0, and LHC ^ 0. After this call, which 
sets the internal variable ^><^=0, if one sets EP (or EP(1))=0 
one has complete control over the integration, and no changes will 
be made in KQ or H unless the user makes them, (The reason HMAX 
must be set to zero, is to protect those users who set EP=0 by 
accident.) If this option is used, one should not change KQ(I) 
by more than ±1 on any step. One should not use this option v?hen 
the integration is being started. 

Anytime the user takes over the control of H, we suggest that 
he get the output provided with L00 of subsection III D on a 
sample problem, as a partial check on the validity of what he is 
doing. 


H, Locating Zeros of Arbitrary Functions (GSTOPs) 

A GSTOP is an indication to the user that a specified function, 
G(J), of the variables T and Y has a zero. The option described 
here“is' used to locate GSTOPs. Applications include the location 
of special output points, and the detection of singularities in 
F which must be found before an attempt to compute F is made. 

Five new parameters are required for this option. 

EG is the number of functions G(j) to be checked for 

zeros. If DAGSl is called with NG ^ 0, this is considered 
a fatal error, and IFLAG is set equal to 8. 

EST0P indicates the index of the function G with the zero 
when a GSTOP is found. Thus G(NST0P) >=« 0. 

LTGS is a vector which indicates how the GSTOP is to be 

located, 

LTGS(J)=0 means do not test for' a zero in G(J), When equal to 
zero, LTGS(J) can be changed at any time, and LTGS(j) 
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can be set equal to zero at any time. 

>0 means locate tiie zero using interpolation, 

<0 means locate the zero using extrapolation. (Note that 
interpolation requires only half as many evaluations 
of G(J) as extrapolation, and gives much better control 
on the error. Thus LTGS(j) < 0 should only be used if 
extrapolation is a necessity.) 

G is a vector of the functions whose zeros' are to be 

located. The user supplies the code to compute G given 
T and Y, 

GT is a vector used for temporary storage of past values 

of G. 

To use this option, the following changes should be made in 
the sample program of subsection III A, 

One should add the following to the type and dimension infor- 
mation, where k & NG, and everything can be scalar if KG=1. 

INTEGER NG,NST0P,LTGS(k) 

D0UBLE PRECISION G(k),GT(k) 

At A— the values of NG and LTGS should be assigned. 

At A^^ if one wants to test for GSTOPs from the beginning, or at 
any other place if one wants to start testing for GSTOPs at a 
later time, insert 

CALL DAGS ( T , Y , IFLAG , NG , NST0P, LTGS , G , GT ) 

This call serves to set up communication. It should only be made 
again if the value of NG is changed. 

The statement at A_ is replaced with (The G. are defined 
below. ) 

A^ G0 T0 IFLAG 

where X^=X 2 =Gij if LTGS(J) < 0 for any J, 

and X^=I^, X^=l^ if LTGS(J) s 0 for all J; and 
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if LTGS(j) > 0 for any J, and otherwise need 

only he included in the list of statement numbers at 

if DV0A0 (see subsection III D) is called with LSS > 

0, in which case X_=I . (Note that one must not call 
o y 

DV0A0 with LSS < 0 if LTGS(j) > 0 for some J. Also, 
if one has called DV0A0' with LEX >0, then one should 
replace the in statement A^ above, with Gg.) 


Between statements labeled A and I, insert the following. 

(It is understood that if LTGS(j) ^ 0 for aJl J, then the statement 
at can be anything and similarly if LTGS a 0 for all J, the 
statement at Gj^ can be anything. ) 



^6 


Compute G(J) for all J which have LTGS(j) > 0 
CALL DAGSl 

G0 T0 ( Ij^, Ig, I^, Ij^,G^,G|j^,G^, l0, I^,Gg) , IFLAG 
Compute G(J) for all J which have LTGS(j) < 0 

G0 T0 G^ 

Do anything special you want to do at the end of each step, 

G0 T0 Ag 

The GSTOP has been found, G(NST0P)f«O. Output results 
and/or whatever else you want to do. If one wants to 
change the definition of G at this time, one can do so 
by storing the values obtained with the new definition 
into the appropriate locations of GT. This is the only 
time the contents of GT should be altered. 

G0 T0 Ag or ST0P or go someplace to restart the integration, 
etc. 

G(NST0P) changes sign, but DAGS could not get good 
convergence when iterating for the zero. This probably 
is due to a bug. 

Print message indicating problem and 
ST0P or G0 T0 G^ 
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I. Extra Output for Checking Out a Program 


When checking out a program it is frequently useful to see 
the values of all the variables involved in the integration. To 
s.et up the subroutine which provides labeled output of the inte- 
gration variables, simply 

CALL DCHK(liffiQ,T,Y,E,KD,EP, IFLAG,H,HMIItA,DELT,TFINAL,MXSTEP, 
KSTEP,KEMAX,EMAX,KQ,YW,DT, I0-) 

where all variables are the same as those in the call to DV0A, 
except for the additional variable 10 which 'is an INTEGER 
variable with usage as follows. 

10 =1 Print everything in the. calling sequence 

=0 Do not print anything 

=-l Print everything in the calling sequence except DT 

If |l0j > 1, then 10 is treated as a vector with dimension at 
least as large as NEQ, and is used as follows. 

I0(l)^ print all information associated with the first equation 
I0(l)=-2 same as 10(1)2:2, exc.ept DT(K,1), k=l,2,... is not 
printed 

I0(l)s-3 do not print variables associated with the first equation, 
for I > 1 when jl0(l)j > 1 , 

I0(I)>O print all information associated with the I-th equation 
I0(l)=O do not print variables associated with the I-th equation 
I0(I)<O same as 10(1) > 0, except DT(K,I), K=1,2,..., is not 
printed 

If one has variational equations and wants output for these 
also, then NTE should be substituted for NEQ in the call, where 
NTE is the total number of equations, and if 10 is a vector it 
should have dimension at least as large as NTE. 

The output is obtained at any time by the simple statement 
CALL DCHKl 

In most cases when this feature is used, one will also want 
to get the output available through the entry DV0A0 described in 
subsection III D. 
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J. Differential Equations With Side Conditions 

In some applications, one is interested in solving eq. (l.l) 
with a side condition of the form 


5(t,y) = 0 


(3J.1) 


where the true solution of eq. (l.l) satisfies ( 3 J.I), and the 
vector cp has dimension the dimension of y. The side condition 
( 3 J.I) is imposed to insure that the numerical solution of eq, (l.l) 
does not drift too far away from satisfying ( 3 J.I). Such" 
problems can be solved using the capability described in subsection 
III E above. Set IfVE(l) = -1, and WVE(3) = 0, where the 0 indicates 
that there are no extra equations to be introduced and that the 
user would like a special return with IFLAG=10 at the time a 
return would ordinarily be made if extra equations were being 
introduced. 

When IFIiAG=10, the user is free to modify y as he wishes in 

order to satisfy (3J.1). In many cases a good solution to (3J.1) 

^ ^ ^ ^ ^ 

will be obtained by setting y = y +6y, where y is the value of y 

c c 

return ed b y the integrator when IFIAG=10, fiy is the vector of 
minimum length which satisfies 


(3^/9y)6y = -^(t,y^) 


(3J.2) 


and (B^/dy) is the Jacobian matrix of ^ with respect to y. 
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IV. USAGE IN THE CASE OF STIFF EQUATIONS 


In tills section the stiff integrator is described following a 
format similar to that in Section III. To simplify the exposition we 
introduce 

u = components of z in eq. (1.1) for which jjsd in eqs. (1.4)-(l,8) 

V = components of z in eq. (l.l) for which A=d, but which may have 

< d at some time in computation. 

w = components of z in eq. (1.1) for which < d ' ' 

X = components of z in eq. (l.l) for which it is possible to have 

A < d at some time in the computation. 

“4 ^4 *4 *4 “4 * 

The corresponding components of y and f are denoted by y ,y ^y ,y , 

uvwx 

f ,f ,f , and f . In all cases the variables which make up the 
u’ v’ w’ X ^ 

individual components of these vectors are in the same order as they 

-» -+ 

are in z, y, and f, respectively. In the usage below, the computations 

of f^ and f^ are placed together, since this simplifies the logic for 

the user. Better results would be obtained by only computing f at 

-4 _> w 

these places, and computing f at those places f is computed. It is 

-♦V -+ -4 u 

also permitted to compute f where f and f are computed and to 

-4U vw ^ 

skip the computation of f^ at the places where the evaluation of f^ 
is indicated; although, of course, one will not get as good results 
doing this, 

A. The Basic Framework 

The stiff integrator requires all of the parameters listed in 
Subsection III A together with those listed below. Usage of the 
parameters given in III A is the same as indicated there, except 
for certain additions noted below. 

IFLAG After DV0A1 is called, usage is the same as in III A 

except when IFLAG=1 or 2 in which case one proceeds as 
follows , 
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=1 Y has been predicted. Compute and (see eq. (4,1)), 

and then call DSTFl. 

*+ 

=2 y has been corrected. Compute and and call DSTFl. 

After DSTFl is called, IFLAG indicates what is to be done, as 
follows 

=1 The first correction of y has been made. Compute f 

W U 

and call DV0A1, 

-4 

=2 The second correction of y has been made. Compute f 
and call DV0A1. 

=3 An extra correction of y is required. Compute f and 

f and call DSTFl. 
w 

=4 Compute partial derivatives, FPS = 5f /3y , and call 

W Vr 

DSTFl. 

~3 A change in the value of JL, see eqs. (1.4)- (1.8), is 
planned. The user may, if he wishes, override any of 
these changes, or he can introduce additional changes 
in jJ, (in many cases, when a change in A is indicated, 
the user may want to insure that the value of is 
"hanged for all equations of a certain class.) The 
procedure for the user to indicate what he wants here 
has not yet been completely specified. Finally call 
DSTFl. 

=6 IFIiAG=6 will only occur if JST0R=O (see subsection IV E 
below), or this type of return is requested through the 
entry DSTF0, (see subsection IV C below). The matrix 
used in the iterative solution of eqs. (1,6)- (1.8) is 
to be changed when IFLAG=6. (This matrix is always 
changed after IFLAG=4, but no return with IFLAG=6 is made 
in this case., since whatever the user wants to do when 
IFIiAG= 6, he can do after computing the Jacobian matrix 
when ZFIiAG= 4.) If JST0R=O, the user may want to do some 
sort of set-up whenever this matrix is changed. Reasons 
for wanting this return when JST0R/O are considered in 
subsection IV C below. 


36 


JPLr Technical Memorandum 33-479 



=7 Possible only when JST0R=O. The matrix formed when 
IFLAG=4 or 6 is to be used to get a correction to 

w 

given the residuals in eg. (1.8). (More on this in 
subsection IV E below. ) Then call DSTFl. 

=8 Fatal error. The value of KEMAX gives the source of 
the error as indicated in subsection III A. If either 
DV0A1 or DSTFl is called, an error message which in- 
cludes the value of KEMAX is printed, and the program 
stopped. 

=9 It appears that the Jacobian matrix is being computed 
incorrectly. To continue: CALL DSTFl, (We recommend 

a careful look at the Jacobian matrix before proceeding.) 

KQ is used as in Subsection III A when in eqs. (l.^)- 

(1.8) . When A < d, KQ, is less than zero, the value of 

g for the I-th equation is given by M0D(-KQ(I),1OO), and 
the value of d-f for the I-th equation is given by 
1-(KQ(I)/100). 

KS vector used to indicate whether the I-th equation is 

implicit, and to indicate the value of A in eqs. (1.4)~ 

( 1 . 8 ) . 

KS(l)<0 means the I-th equation is implicit, otherwise the I-th 
equation is explicit, 

1ks(I) 1<10 means the initial value of KS(l) is never to be changed. 

M0D( |KS(I) j ,10) gives the value of d-A for the I-th equation. 

(Note that one must have d-A ^ and if KS(l) < 0, 
then A < d. If the initial value of KS(l) does 
not satisfy these conditions, IFLAG is set =8.) 
Thus, for example , 

KS(l) =0 the I-th equation is never stiff (i.e, A^). 

=1 on the I-th equation A=d-1 and is not to be changed, 

=10 when starting, the I-th equation is not stiff (A=d), but 
this can change as the integration proceeds. This is 
frequently the best value of KS(l) when starting an 
integration. 
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=-l the I-th equation is implicit, j&=d“l and is not to be 
changed. 

The initial value of KB(I) should satisfy [ks(I)1 < 20. When 
jKS(l)] s 10, the value of jKS(l)| is incremented by multiples 
of 10 in the process of selecting a new value for jJ. 

ASTIF is a logical variable which indicates if any equation 
is stiff ( jJ < d) 

= .TRUE. A < d in some equation 
= .FALSE. == d in all equations 

The user need not set (and shouldn't change) the value 
of ASTIF. 

JST0R is an integer used to indicate how the Jacobian matrix 
is to be stored in FPS. 

=0 FPS can be anything, since it isn't used by DSTF. 

Special returns to the user are made whenever information 
about the Jacobian would ordinarily be required. See 
subsection IV E below. 

=1 FPS(l,J) contains Bf(l)/?y(j), where f(l) = f^ and 
'V'(J)- yj 

=2 FPS(l,j) contains 9f (l)/i5y (J), where f (l) is the I-th 
component of f and y (j) is the J-th component of y ; 

XX X 

see eq, (4,1). JST0R=2 is frequently best when some of 
the KS(I)=0, and the Jacobian is computed analytically. 

S3 FPS(l,j) contains ^f^(l)/^y^( J) ; see eq, (4,1), 

JST0R s 3 is especially useful when one does not know 
which equations are going to be stiff and the Jacobian 
matrix is computed numerically, 

<0 indicates that the Jacobian matrix is a band matrix, 

JST0R is a vector of dimension 3j where JST0R(l) 
indicates the way partials are stored in FPS, JST0R(2) 
gives the number of elements to the left of the main 
diagonal (defined below) and JST0R(3) gives the number 
of elements to the right of the main diagonal. 


38 


JPL Technical Memorandum 33-479 



The main diagonal of the Jacobian matrix consists of 
partials of the form: 

Sf(l)/ay(l), I such that y(l)=z(l), if jjST0R(l)!=l, 

3f^(l)/^^(l), I such that y^(l)3c(l), if |jST0R(l) |=2, 

?f^/I)/Sy^(I), I such that y^(I)^^(I), if 1 jST 0R(1) j^3. 

A 

In the case of first order equations 1=1, In the des- 
cription below we use the abbreviated notation, J = 

Jj 

JST0R(2), Jg = JST0R(3), J^ = + 1 (= the band 

width), and n, n^, n^ are the number of components in 
y> y > y . one must have JST0R(2) and JST0R(3) s 0, 

X MT 

and JST0R(2) + JST0R(3) < n, n^, or n^ depending on the 
value of JST0R(1). 

JST0R(1)=-1 FPS(I,J) contains Bf(l)/Sy(K), where 


( 


J 


K = 


J J+I-Jj^-1 

J+n-J„ 

a 


A 


if 

L 

A 

if J^<Isn-J 

Jj K 

A 
if 

K 


J— 1}2, , , . J Jgj 
J — ’1,2, , ■ • , J^, 
J=lj2, , , , > 


where I is defined as above. Equivalently, one can 
say that Sf(l)/Sy(K) is contained in FPS(I,J), where 


J = 


K 

K-I+Jj^+1 




K^Jg-n 


A 

if i<r^ 

jj 

A 

if J_<l5n-J 

Li K 

A 

if IC>n-J 

K 


K-1,2, , , , J Jg3 

A A ^ 

K— I- , , , , I+Jj^, 

K=n— Jg*^l, • , 4 ^n. 


JST0R(l)=-2 
JST0R( 1)^-3 


Same as liST0R(i)=-1, except f, y, and n are replaced 


with f 


x’ ^x* 


and n 


X 


Same as JST0R(1)=-1, except f, y, and n are replaced 

with f , y , and n . 
v’ ■'w^ w 

Note that the above formulas do permit more than J„ 

A ^ 

elements to the right of the main diagonal when 
and more than elements to the left of it when 
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FPS is a two dimensional array used to store partial deriva- 

tives. See usage of JST0R above, for the way that 
derivatives are stored in FPS. 

IFPS gives the maximum ntimber of rows in FPS, i.e. FPS is 

dimensioned FPS(lFPS,?). IFPS must be at least as large 

■+ 

as the number of components in f, f^, or f^ depending 
on the value’ of JST0E. The number of columns in FPS must be 
s the niimber of components in y, y^, or y^ when JST0R = 

1, 2, or 3, and s Jg = JST0R(2)+JST0R(3)+1 when JST0R(1) < 

0. When JST0R = 0, FPS and hence IFPS can be anything. 

WS is a vector of dimension IWS (see below) used for 

working storage. See subsection IV E for information on 
what is stored in WS. 

OTS is the dimension declared for WS. One must have s 

n +k., where k .=m if JST0R=O, k.=m *(m +3) if JSTOR > 0, 

and k.=2*m *(j +J +2) if JST0R(1) < 0; n = the largest 
J w 3 ^ w 

dimension assumed by y : m = the largest dimension 

w w 

assumed by f^; = JST0R(2), = JST0R(3) if KD=1, and 

with higher order differential equations and give 

the number of equations associated with variables to 

the left and right of the main diagonal. The precise 

definitions of Jg and are somewhat complicated by the 

fact that, as in the last sentence in the usage of 

JST0R, more than variables are permitted to the right 

of the main diagonal when I is close to 1 and more than 

Jg are permitted to the left of the main diagonal when 

I is close to m^. Thus, = the maximum value of K (for 

any l) for which c)f^(l)/Sw(l+K+inax{ 0 ,J 2 -I+l}) is contained 

in FPS and J 2 = the maximum of K for which 

\l-K-max(0,J,-n +i}) is contained in FPS, 
w 

where d is the order of the differential equation associated 
with w(l-K-max( 0 ,J 2 ~n^+l]). 

IW A vector or working storage with dimension s 2*m , where 

m = the largest dimension assumed by f . IW is not 
w -y 

used and thus can be anything if JST0R=O. 
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The sample program for use of the stiff integrator is given 
below. 


In the -type and dimension information given below, (?) indicate 
that the dimension depends on how the parameter is used (see just 
above for parameters which are used only for stiff equations and 
subsection III A for the other parameters), m ^ KEQ, and n & total 
order of the system. 


INTEGER 

LAICAL 

REAL 

DOUBLE PRECISI^i&T 
D0UBLE PRECISION 


NEQ,KD( ? ) , IFLAGjMXSTEP, I^TEP,KEMAX,KQ(m) ,KS(m) , 
JST0R( ? ) , IFPS , IWS , IW( ? ) 

ASTIF 

EP(?),HMIM,EMAX 

F(m) ,DT(20,m) ,FPS(IFPS,? ) ,WS(lWS) 

T,Y(n) ,H,DELT,TFINAL,yN(n) 


Aq Assign the values for lilEQ,KD,MXSTEP,KS, JST0R, IFPS, BfS, 

EP,HMIM,DELT, and TFIUAL 

Assign initial values for T,Y, and H (These will be 
updated by the integrator. Also, depending on the 
initial values in KS, KS may be changed. ) 

Values need not be assigned for IFLAG, ESTEP, KEMAX,KQ, 
IW, AST IF , EMAX, F , DT , FPS , WS , YN 

A^ If certain special options are used, they should be 

set up at this point, 

CALL DSTF(NEQ,T,Y,F,KD,EP, IFLAG ,H , HMINA, BELT , TFINAL , 
MXSTEP, ESTEP, KEMAX, EMAX, KQ,, YN, DT , KS , ASTIF , JST0R , FPS , 
IFPS,WS,IWS,IW) 

G0 T0 A^ 

A^ CALL DV0A1 

A^ G0 T0 (I^,l2,...,lg), IFLAG 

Ij^ C0NTINtlE (special code is required here if certain 

options are used.) 

Ig Compute the f and f parts of F (see eq. (4.1)) 

IF (,N0T. ASTIF) G0 T0 S^ (This statement is optional.) 
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CALL DSTFl 

G0 T0 (3^,82 , IFLAG 
C0NTIWUE 

Compute the part of F (see eq.. (4,1)) 

G0 T0 

Compute partial derivatives, FPS (Only 5f^/dy^ is required, 
but one is permitted to compute more; for example, 

G0 T0 Sq 

The integrator is planning to change see eqs, (l,4)- 
(1,8). At this point the user can do nothing, or 
specify other changes in jj, or override the changes 
planned by the integrator. For instructions, see usage 
for IFLAG=5 above, 

G0 T0 Sq 

If JST0R=O, see subsection IV E below for what to do. 

If DV0A0 WAS CALLED WITH LMF > 0 (see subsection IV C 
below), do whatever you wish. If neither of these 
options is used, then a return with IFLAG=6 should not 
occur. 

If JST0R=O, see subsection IV E below for what to do. 

If JST0R0O, then a return with IFLAG=7 should not occur. 

Print an error message indicating that Jacobian is not 
being computed correctly, and 
ST0P or G0 T0 


Same as in subsection III A. 


END 
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B. Miscellaneous Features Used As in Section III 


Tiie relative error test can be used exactly as described in 
III B. 

To restart an integration, everjrthing in III C applies, except 
that DSTF is substituted for DV0A, and DSTFR is substituted for 
DV0AR. 

Tfie special returns and optional output are available 
exactly as described in III D. 

Tne direct interpolation works as described in III F. Hote 
that anytime it is permitted to CALL DSTFl and anytime a return 
is made from DSTFl, if DAINTl is called, then T, Y, and F must be 
saved and then restored if one intends to continue the integration. 
One change is req,uired when storing information for later 
retrieval due to the fact that KQ(l) may be < 0, In statement 
labeled S^ on page 27 substitute the following for KQ(I); if 
a s d-1 on all equations, substitute ABS(KQ(I))j if jJ may be 
< d-1, then substitute ABS(M0D(KQ(I)),1OO)). 

One gets more direct control over the integration exactly 
as described in III G. 

The GSTOP feature is inserted into the sample program which 
uses DSTF, exactly as is described in III H. 

The extra output described in III I is available exactly as 
described there. 

Differential equations with side conditions are treated much 

as described in III J, except that one uses the variational 

equation entry described in IV D below, and y is modified when 

IFLAG=10 after the call to DSTFl. If this is done, one should 

compute f at the same time as f instead of with f , if it is 
V u w’ 

possible to do so. 
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C. Optional Output and Special Returns for Stiff Equations 
At any place in the program one may insert 

CALL DSTFiZ!(LMF,L0l) 

The variables LMT and L0I are both type INTEGER, and both refer 
to special options as described below. If the user changes one 
of these parameters, DSTF0 must be called for the change to take 
effect, 

LMF determines if a special return is made with IFLAG=6 

whenever the matrix used in the iteration to solve eqs, 
(1.6)~(1,8) is going to be factored. (This return is 
not made if the factorization follows Immediately after 
a Jacobian evaluation. ) This return enables one to 
keep track of the number of matrix factorizations, gives 
one a chance to compute the Jacobian before every 
factorization (a good idea if the Jacobian is extremely 
simple to compute) and makes it possible to use the 
same storage for, EPS and WS, see the last paragraph in 
subsection IV E below. 

>0 Return with IFLAG=6 when the matrix is to be factored 
without an immediately preceding Jacobian evaluation. 

=0 Wo change, (The nominal internal flag is set for no 
return. ) 

<0 Wo return with IFLAG=6 when the matrix is to be factored. 

L0I is used to specify that results connected with the 

iteration to solve eq. (l.6)-(1.8) are to be printed. 

This information is sometimes useful when debugging a 
program. The output to be provided is not yet com- 
pletely specified. 

>0 Give the output. (The size of L0I may determine the 
amount of output. ) 

=0 Wo change (The nominal internal flag is set for no 
output . ) 

<0 Wo output . 
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D, Variational Equations (and Other With Similar Characteristics) 

As in subsection III E, it is convenient to consider the case 
of one additional set of differential equations first. When there 
are stiff equations it is desirable to have the user label whether 
this additional set of equations is a system of variational equations 
as in eq. (3E.2) or is some more general case as permitted by eq. 
(3E.1). In the former case, the same Jacobian matrix is used for 
the variational equations as for the original set of equations, 
while in the latter case the Jacobian matrix 9f /By„ must be stored, 
where y^ is defined as in eq. (3E.1) and f^ = (^’jn+l’ * * * * 

The parameters I3VE and EPSY are defined as in subsection III E. 

In addition, the following parameter is required. 

IREF is used to indicate whether another Jacobian is required. 

=1 The extra equations are variational equations and thus 
a new Jacobian is not required. In this case one must 
have liTVE = i*NEQ where i is a positive integer. Also, 
one should start with KS(l+j*MEQ) = KS(I), 1=1,2, ... ,MEQ, 
j=l,2,. ., ,NVE/NE<^. This is to insure that the same 
value of f, see eqs. (l.4)-(l.8), will be used on all 
equations which utilize the same row of the Jacobian 
matrix. The integrator will not change the value of i 
used on a given equation without changing the value of 
i for all associated variational equations, and vice 
versa. The user who initiates changes in Aj should do 
the same. 

=0 None of the additional equations are permitted to be 
stiff (i.e. have A < d). 

<0 The extra equations are not variational equations, and 
-IREF gives the index in FPS where the first element of 
the Jacobian matrix for the additional equations is 
stored. Thus, for example, if one is using m rows and 
n columns of FPS for the Jacobian matrix of the first set 
of equations and wants to store the Jacobian matrix for 
the extra set of equations in the next available locations. 
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then IREF = -(mn+l). In addition to the extra storage 

required in FPS when ISEF < 0, one must also provide 

extra storage in JST0R and IFPS, both of which are now 

vectors, and in WS and IW. Let ;i]_ last location 

of JST0R used for the first set of equations. Thus if 

JST0R(1) s 0, and if JST0R(1) < 0, dj^=3. Then the 

value of JST0R(j^+l) indicates the way that information 

is stored for the extra set of equations in exactly the 

same way that JST0R(1) indicates the way information is 

stored for the first set of equations. Of course, if 

JST0R(jj^+l) < 0, then JST0 R(Jj_+ 2) and JST0R(J^+3) are 

used just as JST0R(2) and JST0R(3) are used when 

JST0R(1) < 0. IFPS(2.) gives the maximum number of rows 

in the Jacobian matrix for the extra equations in exactly 

the same way that IFPS(l) gives the maximum number of 

rows for the Jacobian matrix of the first set of equations. 

Thus the partial derivatives in the first row of the 

Jacobian for the extra set of equations are stored in 

FPS(-IREF), FPS(-IREF+IFPS(2)), FPS(-IREF+2*IFPS(2) ) , 

, , , . To indicate the extra storage required in WS and 

IW, let m = the largest number of equations in the 
w ^ 

added set which may be stiff at a given time, and n^ = 

the largest number if dependent variables (from y ) which 

may correspond with a stiff equation. The extra storage 

A ^ 

required is just what one might expect: for WS, n 

/\ /\ ” ^ 
where k.=m if JST0R(jT+l) =0, k. = m *(m +3) if 

JST0R(j^+l)>O, and kj=2*m^*( J 2 +J 2 + 2 ) if JST0R( j^+l)<0, 

where J^ and J^ are defined as in the usage of IWS in 

subsection IV A except with JST0R(jj^+2) and JST0R(j^+3) 

substituted for JST0R(2) and JST0R(3)j and for IW, 

2*m . 
w 

Of course, the vectors Y and F are defined as in III E. To 
use this option the following changes should be made in the sample 
program given in subsection IV A. (Changes are similar to those 
given in subsection III E. ) 
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In the dimension information one must have m ^ NEQ+IiVE, and 
n k the total order of the system including the equations which 
have been added. With mixed orders KD must have dimension > m, 
otherwise KD can be a scalar. One should add N7E and lEEF to the 
variables specified to be of type INTEGER, and EPSV to those 
specified to be of type REAL. 

At Aq initial values for WE, EPSV, and IREF must be assigned 
At A^ (before the call to DSTF) 

CALL DSTFV(WVE, EPSV, IREF) 

Code from A^ to the G0 T0 A^ Just below S2 should be replaced 
with 

G0 T0 IFLAG 

I^ C0WTIMUE (Special code is required here if certain 

options are used- ) 

I^Q Compute the f^ and parts of F(l), for 1 i I ^ MEQ 

Sq call DSTFl 

G0 T0 (S^,S2,I^0,Sij,S^,Sg,S.^,Ig,S^,S^Q,S^^), IFLAG 
S^ C03WINUE 

Compute the f^ part of F(I), for 1 ^ I < MEQ 
G0 T0 A^ 

Sj^Q Compute subexpressions which depend only on y, see eq. 

(3E.1) (o?(t,y) and b(t,y) in the case of variational 
equations ) 

Ig Compute f^ and f^^ parts of F(i), I > NEQ 

G0 T0 Sq 

S^3_ C0KTINUE 

Sg Compute part of F(l), I > NEQ, 

If (lFLAG-2) any, S^, Ag 

Note that if F(l) depends only on y (thus not on y ) for I > NEQ, 

3 . 

then none of the equations should be treated as stiff (so have 
KS(l)sO for I > NEQ and IREP=0), and the calculation of the f^ part 
of F(l) at S^can be skipped when IFLAG=2. 
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If one is not using the end of step return (LSS) in subsection 
D, simply set 1 ^^ 2 * 

The more general case in eq.. (3E.4) requires that IREF be a 
vector, and that NVE and EPSV be defined as they are just below 
(3E.4). The K-th additional set of equations is treated as 
indicated by IREF(K). 

IREF(K)=j>0 indicates that the K-th extra set of equations, is 
a set of variational equations associated with the 
j-th set of equations. (The (j+l)-st set of equations 
= j-th extra set of equations.) In this case one 
must have ]WE(K+2)=i*tj,^ where i is a positive 
integer and n =NVE(j+l) if j>l and =NEQ if j=l. 

As before, the initial values of KS(l) should be the 
same for all equations which use the same row of a 
given Jacobian matrix, 

IREF(K)=0 indicates that there are no stiff equations in the 
K-th extra set of equations. 

IREF(K)<0 indicates that the K-th extra set of equations is 
not variational equations, and specifies the 
starting point in FPS where the Jacobian matrix is 
to be found as described under IREF<0 above. As 
before, the way the Jacobian is stored in FPS is 
indicated in JST0R starting with JST0R( jj^+1) , where 
jj^ is the last value of JST0R used by the K-th set 
of equations. (IREF(K)sO does not use any storage 
in JST0R. ) Also, IFPS(K+l) gives the maximum row 
dimension of the Jacobian stored in FPS. The extra 
storage required in WS and IW is computed using the 
obvious generalization of the formulas given earlier 
under IREF<0. 

IREF must satisfy the following: If IREF(K)<0 then IREF(i)^0 

for all i<K; if IREF(K)=j>0, then j>IREF(i) for all ±<K and 
either j=l, or IREF(j-l)<0. 
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In order to handle several additional sets of equations, one 
should modify the changes given earlier in this subsection by 
dimensioning IREF, and then changing the code from on, as follows. 

C0KTIKUE (Special code is required here if certain 
options are used, ) 

K=1 

G0 T0 I2 

IF (K.EQ.l) G0 T0 I^ 

' IF (lEEF(K-l).EQ.O) G0 T0 

I^ Compute f^ and f^ parts for the K-th set of F(l)'s 

Sq call DSTFl 

G0 T0 (S^,S2,E^,Si^,S^,Sg,S^,Ig,S^,Sj_Q,S^^), IFIAG 

S^ C0NTINUE 

G0 T0 

K=NVE(2)+1 

Compute subexpressions for the K-th set of equations 
(NVE(2),-th additional set) which depend only on ccm- 
ponents of Y from previous sets. 

G0 T0 Ej^ 

“4 

S-,, Compute f part for the K-th set of F(l)’s, 

11 u 

G0 T0 Ag 

Sg Compute f^ part of all f(I)’s (Those components of F 

which depend on]^ on components of Y from an earlier 
set of equations need not be computed at this time. ) 

G0 T0 Ag 

E, Jacobian Has Special Structure (e.g. Sparse) 

By setting JST0R=O, the user can take control over the 

solution of eqs. (l,6)-(l,8). In this subsection, the following 

notation is used, where y , z , and f are defined at the- beginning 
^ •'w’ w’ w 

of Section IV. 
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-y,’ 

c = (Ci>C2>-">V ' 

-+ T "* 

9 " ( 92 ’ 92 ’ * • * > 9 ( 0 ^) ~ 


i=l,2,...,m is ttie value of 1 (see eqs. (l.4)-(l,8)) 

(a) 


used with 




V4e.1) 


is the order of the differential equation corresponding 
to cp^ , 

is the integration order used on the equation corres- 
ponding to 


p. . = 3cD./aT| i=l,2, 
10 10 




y 


The option described here presumably is of interest only if there 

is something special about the matrix (p . . ) which the user wants 
10 

to take advantage of in obtaining the solution of eqs. (1.6)-(l,8). 
An example of such a situation is the case when most of the p 

10 

are zero. We restrict our attention here to the case when the 
user linearizes cp, and solves the resulting linear problem as an 
approximation to the solution of eqs. (l.6)-(l.8). There are 
other ways the user might approach the problem of getting an 
approximate solution to eqs. (l.6)-(l.8). 

Consider first the case when d.=l, and thus n=m, and 
The user solves for a correction, 6?, to E using the formula 

A6§ = r = (^E.2) 
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where A=(a^^), if Vj and a^^=p_-WS(i) , WS(i)=(l/h)g^ 

if the i-th equation is explicit and =0 if the i-th equation is 
implicit, r^=WS(n+i)= the residual that results from subtracting 
the left side of eq, (1,8) from the right side, and the contents 
of WS are computed by DSTF, When IFLA.G=6 and after computing the 
Jacobian when IFLAG=4, the user should compute the matrix A, and 
do any preliminary set up that will be used for later iterations. 
When H’LAG=7, the user should solve eq, (4 e, 2), and store the 
vector 65 in the storage previously occupied by r. In all of 
these cases, one continues by calling DSTFl, 

When some of the d. are greater than one, n=d,+d^+** *+d > m, 

1 12m’ 

and an nxn linear system results from the linearizing eqs. (1,6)- 
(1.8), However, the special structure of eqs, (1.6) and (1.7) makes 
it possible to eliminate n-m variables in eq. (1.8), One is then 
left with the mxm system (4 e, 2) to solve, and DSTF solves for the 
remaining varaibles. Usage is the same as described above, ex- 
cept that the matrix A is more complicated to obtain, and the 
contents of WS require a more general definition. Let 4:^=0, k^=d^ 
+d 2 +»**+d^, ±>0, WS(j) contains a coefficient used in obtaining 
the value of T]^. More specifically, WS(k^ WS(k^) contain 

coefficients used for . , . . . ^ \ The first cells 

starting j^ith WS(k^_^+l) (if jj=0, then no cells) contain 


h 


^ij ^i 


3 • • • ^ 



the next cell contains 


A.-d. 

h d.-jJ. i-th equation is explicit). 


or 0 (if implicit). 
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the remaining cells (none if contain 


h p . , , , . jh 


-(d -Ji -1) 

J. J. . 


1 1 


As before WS(n+i) contains the i-th residual, r^, to eq^. (1.8). 
From the structure of eqs, (l.6)~(l.8) and the fact that the 
initial estimates are such that the residvials for eqs. (1.6) and 
(1.7) are zero, it follows that A can be defined as follows. 


10 


=[ws(k^_^+l)*p(i,k^_^+l)+. . .+WS(k^_^+ji^)*p(i,k^^^+A^)] 


p(i>k +1) 

J d 


if ij^O 


p(i,k. +jJ.+l) - WS(k. ^+^.+1) if i=j 

d d. J J d- d 


d 


+rws(k._ +je.+2)*p(i,k._ +^.+2)+*“+WB(k.)*p(i,k.) . 

L d**-^ d d~*^ d d d - 


Where either of the sums in brackets may not contain any terms 

(and hence be equal to zero) depending on the values of 1. and 

d 

d.. (Uote that k.=k. ,+d.; and p(i,k. +1), . . . ,p(i,k. ) are used 

d d d"*“^ d d""*^ d 

in forming a . . . ) 

^d 

When JST0Rji^O, DSTF selects the value of SL^ automatically, 
but when JST0R=O the value of is only changed if the user 
initiates the change. 

The first n and the next m locations in WS are always used 
w w 

as described above. When JST^Rj^O, the next locations are 

used as working storage by the linear algebra subroutines. And 
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2 

finally, the next locations (if JST0R>O) or the next m^*(2*o{g-l) 
locations (if JST0R(l)<O, and where of denotes the band width in 
the matrix A) are used to store the matrix A and the factored form 
of A. When there are variational equations this same pattern is 
followed for the additional equations, always starting with the 
first unused location in WS* 

If there are no variational equations, one always recomputes 
FPS whenever A is factored (use LMF>0 in the options described in 
subsection IV C and compute FPS when IFIAG=6), and FPS is not 
used for anything but forming the matrix A in eq, (4 e, 2), then one 
can use part of the storage required for WS to hold FPS. Use an 
EQUIVAIOTCE statement to use the same location for FPS(l,l) and 
WS(k), where k^^+1 and j = maximum value ever attained by n^+3*m^. 


F. Control of the Type of Past Information Stored 

Since we do not as yet have a good automatic algorithm for 
changing A when j0=d (see eqs. (l,U)-(l.8) ) , and even if we did 
there would still be cases when the user might do a better job, a 
feature is included which enables the user to change Z at anytime 
in the computation. To use this feature one executes the following 
statement -at anytime 

CALL DSTFC 

and then continues as one would if this statement were absent. At 
the first opportunity DSTF will return control to the user with 
IFLAG=5. The user can then change the vector KS as described in 
subsection IV A under the usage of IFLAG=5. 


G. User Control of Jacobian Evaluations 

If the optional return with IFLAG=6 does not let the user 
evaluate the Jacobian matrix as often as he wants to, he can 
evaluate the Jacobian more frequently using the feature described 
here. At anytime execute 

CALL pSTFJ 
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and then continue as usual. At the first opportunity DSTF will 
return control to the user with IFLAG=4 at v?hich time the user 
should compute the Jacobian matrix and call DSTPl. Note that it 
is permitted to change FPS without the new Jacobian being used in 
the iteration to solve eqs. (l.6)-(l,8), (.We assume that the same 
storage is not being used for FPS and WS.) For example, one may 
want to do this when integrating variational equations. Of course, 
if this is the case one would not call DSTFJ, 


H. Numerical Evaluation of the Jacobian Matrix 

We plan to provide a subroutine DPAPT which the user can call 
anytime he needs to compute partial derivatives for the Jacobian 
matrix. More work is required before specifications for this sub- 
routine can be given. 
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V. 


CONCLUDING REMAEKS 


We plan to use parts of this document in the subroutine write-up 
to be prepared later. However, when this work is completed it is 
likely that the usage will be slightly different, and the write-up 
of the subroutines will include additional material to aid people who 
use them. In addition, a shorter write-up will be available for those 
with simple problems. We also plan to give some thought to providing 
a program which will only require the user to provide the code for 
the derivative calculations, and which will enable the user to specify 
things such as form of output, error tolerance, initial conditions, etc. 
through special input parameters. Of course, such a program would not 
have near the flexibility of the software described here. 

It would have been nice to include the ability to automatically 
select the type of error test to be used. That is, should the local 
error be kept less than some specified e or less than e times some 
function, and if so, what function? The answer, of course, depends 
on the global effect of local errors, and these effects seem impossible 
to estimate automatically without an inordinate amount of computation. 

In [l], we give suggestions on how to select the type of error test. 
Unfortunately, -users avoid reading subroutine write-ups if they can, 
and they frequently can. The result is that users frequently use a 
type of test that results in an inefficient computation. Is there a 
good solution to this problem? 

For completeness, we mention one more feature which will be added 
if there is sufficient interest. This feature will permit the user 
to group his differential equations, with all equations in a group 
integrated with the same stepsize, but with the possibility that 
different stepsizes will be used for equations in different groups. 

All stepsizes v?ould be some power of two times the smallest stepsize. 
This feature woixld introduce an extra version of DV0A (and DSTF) if 
it requires a significant amount of extra code in DV0A, or if it has 
any more than negligible effect on the integration overhead when the 
feature is not used. All of the features including already have 
satisfied these two constraints. 
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Jacobian Generator (DPAET) 

TiS 

I 

I 



Dashed lines are for optional calls, 

I— B I means A calls B, and B returns control to A. 


FIGURE 1: 


Communication Between the Subroutines and the User 
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